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Abstract 

The prevailing theory for the molecular basis of evolution involves genetic mutations that ultimately generate the heritable phenotypic 
variation on which natural selection acts. However, epigenetic transgenerational inheritance of phenotypic variation may also play an 
important role in evolutionary change. A growing number of studies have demonstrated the presence of epigenetic inheritance in a 
variety of different organisms that can persist for hundreds of generations. The possibility that epigenetic changes can accumulate 
over longer periods of evolutionary time has seldom been tested empirically. This study was designed to compare epigenetic changes 
among several closely related species of Darwin's finches, a well-known example of adaptive radiation. Erythrocyte DNA was 
obtained from five species of sympatric Darwin's finches that vary in phylogenetic relatedness. Genome-wide alterations in genetic 
mutations using copy number variation (CNV) were compared with epigenetic alterations associated with differential DNA meth- 
ylation regions (epimutations). Epimutations were more common than genetic CNV mutations among the five species; furthermore, 
the number of epimutations increased monotonically with phylogenetic distance. Interestingly, the number of genetic CNV muta- 
tions did not consistently increase with phylogenetic distance. The number, chromosomal locations, regional clustering, and lack of 
overlap of epimutations and genetic mutations suggest that epigenetic changes are distinct and that they correlate with the evo- 
lutionary history of Darwin's finches. The potential functional significance of the epimutations was explored by comparing their 
locations on the genome to the location of evolutionarily important genes and cellular pathways in birds. Specific epimutations were 
associated with genes related to the bone morphogenic protein, toll receptor, and melanogenesis signaling pathways. Species- 
specific epimutations were significantly overrepresented in these pathways. As environmental factors are known to result in heritable 
changes in the epigenome, it is possible that epigenetic changes contribute to the molecular basis of the evolution of Darwin'sf inches. 
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Introduction 

Epigenetic change has been postulated to play a role in the 
ecology and evolution of natural populations (Richards et al. 
2010; Holeski et al. 2012; Liebl et al. 2013). Epigenetic 
changes are broadly defined as "molecular processes 
around DNA that regulate genome activity independent of 
DNA sequence and are mitotically stable" (Skinner et al. 
2010). Some epigenetic processes are also meiotically stable 
and are transmitted through the germline (Anway et al. 2005; 
Jirtle and Skinner 2007). These epigenetic mechanisms, 
such as DNA methylation, can become programmed 



(e.g., imprinted) and inherited over generations with potential 
evolutionary impacts. Environmental factors have been shown 
to promote the epigenetic transgenerational inheritance of 
phenotypic variants (Skinner et al. 2010). In recent years, the 
importance of environmental cues in the induction of such 
variation has been widely acknowledged (Bonduriansky 
2012). Thus, like genetic change (Greenspan 2009), epige- 
netic change may also play an important role in evolution 
(Guerrero-Bosagna et al. 2005; Day and Bonduriansky 201 1; 
Geoghegan and Spencer 2012, 2013a, 2013b, 2013c; 
Klironomos et al. 2013). 
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In order for inherited epigenetic changes to play a signifi- 
cant role in microevolution, they must persist for tens of gen- 
erations, or longer (Slatkin 2009). It is conceivable that 
epigenetic changes may also accumulate over longer periods 
of evolutionary time, contributing to processes such as adap- 
tive radiation (Rebollo et al. 2010; Flatscher et al. 2012). This 
hypothesis assumes that epigenetic changes persist over thou- 
sands of generations. An initial step in testing this hypothesis 
would be to compare epigenetic differences among closely 
related species, and whether such changes accumulate over 
short spans of macroevolutionary time. For example, do epi- 
genetic changes accumulate with phylogenetic distance? 
Addressing this question was the primary goal of this study. 

The study was designed to explore the relationship between 
epigenetic changes and the evolutionary history of several spe- 
cies of Darwin's finches in the Galapagos Islands. This group of 
birds has been central to work on a variety of important topics 
in evolutionary biology, including adaptive radiation, character 
displacement, rapid evolution, hybridization between species, 
evolutionary developmental mechanisms, and the effect of 
invasive pathogens and parasites (Grant and Grant 2008; 
Huber et al. 2010; Donohue 2011). The adaptive radiation 
of Darwin's finches over a period of 2-3 Myr resulted in 14 
extant species that fill distinct ecological niches. These species 
show striking variation in body size and the size and shape of 
their beaks (Grant and Grant 2008). Darwin's finches were 
selected for study because they are a well-studied example 
of the evolution of closely related species into different eco- 
logical niches (Grant and Grant 2008; Donohue 201 1). 

Natural selection is a process in which environmental fac- 
tors influence the survival and reproductive success of individ- 
uals bearing different phenotypes. Only selection on 
phenotypic traits with a heritable basis can lead to evolution- 
ary change (Endler 1 986). Observations indicate that epige- 
netic mechanisms have a role in influencing genomic 
variability (Huttley 2004; Ying and Huttley 2011). As epige- 
netic changes are also influenced by environmental factors, 
and can be heritable across generations (Skinner et al. 2010), 
they provide another molecular mechanism that can influence 
evolutionary change. Although Lamarck (1802) proposed that 
environmental factors can influence inheritance directly, his 
mechanism has not been widely recognized as a component 
of modern evolutionary theory (Day and Bonduriansky 201 1). 
Recent work in epigenetics shows that epigenetic changes 
can, in fact, increase the heritable phenotypic variation avail- 
able to natural selection (Holeski et al. 201 2; Liebl et al. 201 3). 
Thus, epigenetics appears to provide a molecular mechanism 
that can increase phenotypic variation on which selection acts 
(Skinner 2011). The integration of genetic and epigenetic 
mechanisms has the potential to significantly expand our un- 
derstanding of the origins of phenotypic variation and how 
environment can influence evolution. 

For example, Crews et al. (2007) investigated the ability of 
an environmental factor (toxicant) to promote the epigenetic 



transgenerational inheritance of alterations in the mate pref- 
erences of rats, with consequences for sexual selection. An F0 
generation gestating female rat was exposed to the agricul- 
tural fungicide vinclozolin transiently. A dramatic alteration in 
the mate preferences of the F3 generation was observed 
(Crews et al. 2007) along with epigenetic alterations 
(termed epimutations) in the germline (sperm) (Guerrero- 
Bosagna et al. 2010). Transgenerational transcriptome 
changes in brain regions correlated with these alterations in 
mate preference behavior were also observed (Skinner et al. 
2008, 2014). Thus, an environmental factor that altered mate 
preference was found to promote a transgenerational alter- 
ation in the sperm epigenome in an imprinted-like manner 
that was inherited for multiple generations (Crews et al. 
2007; Skinner et al. 2010). Studies such as these suggest 
that environmental epigenetics may play a role in evolutionary 
changes through processes, such as sexual selection. 

Recent reviews suggest a pervasive role for epigenetics in 
evolution (Rebollo et al. 2010; Day and Bonduriansky 201 1; 
Kuzawa and Thayer 201 1; Flatscher et al. 2012; Klironomos 
et al. 2013). The primary goal of this study was to test whether 
epigenetic changes accumulate over the long periods of evo- 
lutionary time required for speciation with adaptive radiation. 
Genome wide analyses were used to investigate changes in 
genetic and epigenetic variation among five species of 
Darwin's finches. The measure of genetic variation was copy 
number variation (CNV), which has been shown to provide 
useful and stable genetic markers with potentially more phe- 
notypic functional links than point mutations such as single 
nucleotide polymorphisms (SNPs) (Lupski 2007; Sudmant et al. 
201 3). CNVs involve an increase or decrease in the number of 
copies of a repeat element at a specific genomic location. 
Recently, CNV changes in primates and other species have 
been shown to be very useful genetic measures for comparing 
evolutionary events (Nozawa et al. 2007; Gazave et al. 201 1 ; 
Poptsova et al. 2013). CNV changes are involved in gene du- 
plication and deletion phenomena, as well as repeat element 
phenomenon such as translocation events and can be influ- 
enced by DNA methylation (Skinner et al. 2010; Macia et al. 
201 1; Tang et al. 2012). The measure of epigenetic variation 
used was differential DNA methylation sites, which are known 
to be stable and heritable (Skinner et al. 2010). Comparing 
data for both genetic mutations (i.e., CNV) and epimutations 
(i.e., DNA methylation) allowed the relative magnitudes of 
these sources of variation to be compared across the five spe- 
cies included in the study. 

Materials and Methods 

Finch Field Work and Collection of Blood 

Blood samples were collected from birds captured January- 
April 2009 at El Garrapatero, a lowland arid site on Santa Cruz 
Island, Galapagos Archipelago, Ecuador (Koop et al. 2011). 
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Birds were captured with mist nests and banded with num- 
bered Monel bands to track recaptures. Birds were identified, 
aged, and sexed using size and plumage characteristics. A 
small blood sample (90 pi) from each bird was collected in a 
microcapillary tube through brachial venipuncture. Samples 
were stored on wet ice in the field, then erythrocytes purified 
by centrifugation and cells stored in a -20 °C freezer at a field 
station. Following the field season, samples were placed in a 
-80 °C freezer for longer term storage. All procedures were 
approved by the University of Utah Institutional Animal Care 
and Use Committee (protocol #07-08004) and by the 
Galapagos National Park (PC-04-10: #005441 1). 

DNA Processing 

Erythrocyte DNA was isolated with DNAeasy Blood and Tissue 
Kit (Qiagen, Valencia, CA) and then stored at -80 °C prior to 
analysis. DNA was sonicated following a previously described 
protocol (without protease inhibitors) (Tateno et al. 2000) and 
then purified using a series of washes and centrifugations 
(Ward et al. 1999) from variable number of animals per spe- 
cies analyzed. The same concentrations of DNA from individ- 
ual blood samples were then used to produce pools of DNA 
material. Two DNA pools were produced in total per species, 
each one containing the same amount of DNA from different 
animals. The number of individuals used per pool is shown in 
supplementary table S6, Supplementary Material online. 
These DNA pools were then used for chromosomal genomic 
hybridization (CGH) arrays or chromatin immunoprecipitation 
of methylated DNA fragments (MeDIP). 

CNV Analysis 

The array used for the CNV analysis was a CGH custom design 
by Roche Nimblegen that consisted of a whole-genome tiling 
array of zebra finch (Taeniopygia guttata) with 720,000 
probes per array. The probe size ranged from 50 to 75 mer 
in length with median probe spacing of 1,395 bp. Two 
different comparative (CNV vs. CNV) hybridization experi- 
ments were performed (two subarrays) for each species in 
query (Geospiza fuliginosa [FUL], G. scandens [SCA], 
Camarhynchus parvulus [PAR], and Platyspiza crassirostris 
[CRA]) versus control G. fortis (FOR), with each subarray in- 
cluding hybridizations from DNA pools from these different 
species. Two DNA pools were built for each species (supple- 
mentary table S6, Supplementary Material online). For one 
subarray of each species, DNA samples from the experimental 
groups were labeled with Cy5 and DNA samples from the 
control lineage were labeled with Cy3. For the other subarray 
of each species, a dye swap was performed so that DNA sam- 
ples from the experimental groups were labeled with Cy3 and 
DNA samples from the control lineage were labeled with Cy5. 

For the CNV experiment raw data from the Cy3 and Cy5 
channels were imported into R (R Development Core Team 
2010), checked for quality, and converted to MA values 



(M=Cy5 - Cy3; A = [Cy5 + Cy3]/2). Within array and be- 
tween array normalizations were performed as previously de- 
scribed (Manikkam et al. 2012). Following normalization, the 
average value of each probe was calculated and three differ- 
ent CNV algorithms were used on each of these probes 
including circular binary segmentation from the DNA copy 
(Olshen et al. 2004), CGHseg (Picard et al. 2005) and 
cghFlasso (Tibshirani and Wang 2008). These three algorithms 
were used with the default parameters. The average values 
from the output of these algorithms were obtained. A thresh- 
old of 0.04 as a cutoff was used on the summary (average of 
the log-ratio from the three algorithms) where gains are 
probes above the positive threshold and losses are probes 
below the negative threshold. Consecutive probes (>3) of 
gains and losses were used to identify separate CNV regions. 
A cutoff of three-probe minimum was used and those regions 
were considered a valid CNV. The statistically significant CNVs 
were identified and P values associated with each region pre- 
sented. A cutoff of P< 10~ 5 was used to select the final re- 
gions of gains and losses. 

Differential DNA Methylation Regions Analysis 

MeDIP was performed as previously described (Guerrero- 
Bosagna et al. 2010) as follows: 6pg of genomic DNA was 
subjected to series of three 20-pulse sonications at 20% am- 
plitude and the appropriate fragment size (200-1 ,000 ng) was 
verified through 2% agarose gels; the sonicated genomic 
DNA was resuspended in 350 pi TE buffer and denatured 
for 10min at 95 °C and then immediately placed on ice for 
5min; 100 pi of 5x IP buffer (50 mM Na-phosphate pH 7, 
700 mM NaCI (PBS), 0.25% Triton X-100) was added to the 
sonicated and denatured DNA. An overnight incubation of the 
DNA was performed with 5pg of antibody anti-5- 
methylCytidine monoclonal from Diagenode (Denville, NJ) at 
4°C on a rotating platform. Protein A/G beads from Santa 
Cruz were prewashed on PBS-BSA (bovine serum albumin) 
0.1% and resuspended in 40 pi 1x IP (immunoprecipitation) 
buffer. Beads were then added to the DNA-antibody complex 
and incubated 2 h at4°C on a rotating platform. Beads bound 
to DNA-antibody complex were washed three times with 1 ml 

1 x IP buffer; washes included incubation for 5 min at 4 °C on 
a rotating platform and then centrifugation at 6,000 rpm for 

2 min. Beads DNA-antibody complex were then resuspended 
in 250 pi digestion buffer (50 mM Tris-HCI pH 8, 10mM eth- 
ylenediaminetetraacetic acid, 0.5% SDS (sodium dodecyl sul- 
fate) and 3.5 pi of proteinase K (20 mg/ml) was added to each 
sample and then incubated overnight at 55 °C on a rotating 
platform. DNA purification was performed first with phenol 
and then with chloroform:isoamyl alcohol. Two washes were 
then performed with 70% ethanol, 1 M NaCI, and glycogen. 
MeDIP-selected DNA was then resuspended in 30 pi TE buffer. 

The array used for the differential methylation analysis was 
a DNA-methylated custom array by Roche Nimblegen that 
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Fig. 1. — Number of epimutations and genetic mutations in relation to the phylogenetic relationships of five species of Darwin's finches. Photographs (by 
J.A.H.K. or S.A.K.) show variation in bill size and shape. Numbers on branches are the number of differences (three or more probes; table 1) in epimutations 
(DMR; in red) .and genetic mutations (CNV; in blue) for each of four species, compared with a single reference species FOR (asterisk). The phylogram is based 
on allele length variation at 1 6 polymorphic microsatellite loci (from Petren et al. 1 999). The topology of the tree is similar to that proposed by Lack (1 947) on 
the basis of morphological traits. 



consisted of a whole-genome tiling array of zebra finch 
(Taeniopygia guttata) made of four 2.1 M and one 3x720k 
array with 8,539,570 probes per array. Probe sizes were 50- 
75 mer in length and median probe spacing was 200 bp. Two 
different comparative (MeDIP vs. MeDIP) hybridization exper- 
iments were performed (two subarrays) for each experimental 
species (FUL, SCA, PAR, CRA) versus control FOR, with each 
subarray including hybridizations from MeDIP DNAfrom DNA 
pools from these different species (supplementary table S6, 
Supplementary Material online). For one subarray of each spe- 
cies, MeDIP DNA samples from the experimental groups were 
labeled with Cy5 and MeDIP DNA samples from the control 
lineage were labeled with Cy3. For the other subarray of each 
species, a dye swap was performed so that MeDIP DNA sam- 
ples from the experimental groups were labeled with Cy3 and 
MeDIP DNA samples from the control lineage were labeled 
with Cy5. 



For each comparative hybridization experiment, raw data 
from both the Cy3 and Cy5 channels were imported into R, 
checked for quality, and converted into MA values. The nor- 
malization procedure is as previously described (Guerrero- 
Bosagna et al. 2010). Following normalization each adjacent 
>3 probe set value represents the median intensity difference 
between FUL, SCA, PAR and CRA and control FOR of a 600- 
bp window. Significance was assigned to probe differences 
between experimental species samples and reference FOR 
samples by calculating the median value of the intensity dif- 
ferences as compared with a normal distribution scaled to the 
experimental mean and standard deviation of the normalized 
data. A Z score and P value were computed for each probe 
from that distribution. The statistically significant differential 
DNA methylation regions (DMR) were identified and P values 
associated with each region represented, as previously de- 
scribed (Guerrero-Bosagna et al. 2010). 
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Fig. 2. — Number of epimutations and genetic mutations associated with Darwin's finches. The number of differential DMR epimutations and CNV 
genetic mutations {A). DMR and CNV that differ significantly (P< 1CT 5 ) from the reference species (FOR) are presented for all oligonucleotide probes, 
compared with peaks of three or more adjacent probes. The epimutations with an increase (Up) or decrease (Down) in DNA methylation are indicated. Those 
genetic mutations with an increase (Gain) or decrease (Loss) in CNV are indicated. Venn diagrams for epimutations (6) and genetic mutations (O show 
overlaps between epimutations (DMR) and genetic mutations (CNV) among species. The species and total number of sites compared are listed on the outside 
of each colored elliptical. 



Additional Bioinformatics and Statistics 

The July 2008 assembly of the zebra finch genome (taeGutl , 
WUSTL v3.2.4) produced by the Genome Sequencing Center 
at the Washington University in St Louis (WUSTL) School of 



Medicine was retrieved (WUSTL 2008). A seed file was con- 
structed and a BSgenome package was forged for using the 
Finch DNA sequence in the R code (Herve Pages BSgenome: 
Infrastructure for Biostrings-based genome data packages. R 
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Fig. 3. — Phylogenetic distance is correlated with epigenetic changes, 
but not genetic changes. Branch lengths in figure 1 were used as measures 
of phylogenetic distance. The number of epimutations increased with phy- 
logenetic distance (Spearman Rho=1.0, P< 0.0001). In contrast, the 
number of genetic mutations did not increase with phylogenetic distance 
(Spearman Rho = 0.8, P=0.2). 

package version 1 .24.0). This sequence was used to design 
the custom tiling arrays and to perform the bioinformatics. 

The chromosomal location of CNV and DMR clusters used 
an R-code developed to find chromosomal locations of clus- 
ters (Skinner et al. 201 2). A 2-Mb sliding window with 50,000 
base intervals was used to find the associated CNV and DMR 
in each window. AZ-test statistical analysis with P< 0.05 was 
used on these windows to find the ones with overrepresented 
CNV and DMR were merged together to form clusters. A 
typical cluster region averaged approximately 3 Mb in size. 

The DMR and CNV association with specific zebra finch 
genes and genome locations used the Gene NCBI database 
for zebra finch gene locations and correlated the epimutations 
associated (overlapped) with the genes. The three adjacent 
probes constituted approximately a 200-bp homology 
search. The KEGG (Kyoto Encyclopedia of Genes and 
Genomes) pathway associations were identified as previously 
described (Skinner et al. 201 2). Statistically significant overrep- 
resentation uses a Fisher's exact analysis. 

Spearman Rank correlation coefficients were used to test 
for a relationship between phylogenetic distance and epige- 
netic and genetic changes (Whitlock and Schluter 2009). 

Results 

Phylogenetic relationships of the five finch species in this study 
are shown in figure 1 . The taxa chosen for this study included: 



Two species of ground finches, FOR and FUL, which have 
crushing beaks with relatively deep bases; the cactus finch 
SCA, which has a long thin beak used for probing flowers; 
the small tree finch PAR, which has curved mandibles used for 
applying force at the tips; and the vegetarian finch CRA, 
which has a relatively short stubby bill used for crushing 
food along its entire length (Grant and Grant 2008; 
Donohue 2011; Rands et al. 2013). FOR was selected as a 
reference species for comparing genetic and epigenetic alter- 
ations among the remaining four species. Branch lengths in 
figure 1 were used as measures of phylogenetic distance. 

The experimental design used purified erythrocytes from 
the different species. Although DNA sequences are the 
same for all cell types of an organism, the epigenome is dis- 
tinct for each cell type, providing a molecular mechanism for 
the genome activity and functions that differ among different 
cell types (Skinner et al. 2010). Therefore, to investigate the 
overall epigenome requires a purified cell type. As birds have 
erythrocytes (red blood cells) that contain nuclei, samples of 
purified erythrocytes were collected from each of the Darwin's 
finch species to obtain DNA for molecular analysis. 

The epigenetic alterations termed epimutations were as- 
sessed through the identification of differential DMR. The 
DMR were identified with the use of MeDIP with a methyl 
cytosine antibody, followed by a genome wide tiling array 
(Chip) for an MeDIP-Chip protocol (Guerrero-Bosagna et al. 
2010). Although other epigenetic processes such as histone 
modifications, chromatin structure, and noncoding RNA are 
also important, DNA methylation is the best known epigenetic 
process associated with germline-mediated heritability and en- 
vironmental manipulations (Skinner et al. 2010). Genetic var- 
iation was assessed using CNVs (i.e., amplifications and 
deletions of repeat elements) in the DNA using a CGH proto- 
col (Pinkel and Albertson 2005; Gazave et al. 201 1). 

The reference genome used for the analysis was that of the 
zebra finch (Taeniopygia guttata) (Clayton et al. 2009), which 
had a preliminary estimate of greater than 83% similarity with 
a partial shotgun sequence of a Darwin's finch genome 
(Rands et al. 2013). This study actually suggests a much 
higher degree of identity. The zebra finch genome was tiled 
in a genome wide array with a 200-bp resolution and for a 
CGH array with a 1 ,500-bp resolution. These arrays were used 
in a competitive hybridization protocol between FOR (refer- 
ence species) and the other four species (Guerrero-Bosagna 
et al. 201 0). Differential hybridization using two different fluo- 
rescent DNA labeling tags identified the CNV with CGH using 
genomic DNA and the epimutation DMR with a MeDIP-Chip 
protocol. A statistical significance threshold of P< 10~ 5 was 
set for the CNV or epimutation to be identified as a gain or 
loss compared with the reference species (fig. 2 and supple- 
mentary tables S1 and S2, Supplementary Material online). 
The data for all probes (oligonucleotides on the arrays) are 
presented. However, the criteria used to identify the CNV 
and DMR required the involvement of three or more adjacent 
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Fig. 6. — Chromosomal locations for clusters of CNV and DMR. The chromosome number and size are presented in reference to the zebra finch 
genome. The chromosomal location of statistically significant (P< 1CT 5 ) overrepresented clusters of CNV (A) and DMR (6). The legend shows species and 
total number of clusters. 
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probes in the genome sequence having significant differential 
hybridization. These selection criteria reduce the number of 
false positives and provide a more reliable comparison (fig. 2). 
Therefore, the data presented used stringent criteria and rep- 
resent the most reproducible epimutations and genetic CNV 
mutations among all three different experiments. 

The increases or decreases in DNA methylation for the DMR 
are presented, along with the total number of epimutations in 
figure 2. The majority of epimutations for all species but FUL 
involves a decrease in DNA methylation (fig. 2/4). The gains or 
losses in CNV are also presented, along with the total number 
of genetic alterations. The majority of genetic mutations for all 
species but PAR involves an increase in CNV number. 
Interestingly, the number of epimutations observed was gen- 
erally higher, using the criteria selected, than the number of 
genetic alterations (fig. 2). However, the overall magnitude of 
epigenetic change was comparable to that of genetic change. 
Data for the five different species are shown in figure 1 for 
both epimutations (red) and genetic alterations (blue). The 
number of epimutations was significantly correlated with phy- 
logenetic distance, whereas the number of genetic mutations 
was not (fig. 3). 

The chromosomal locations of the CNV for the different 
finch species are shown in figure 4. CNVs were found on most 
chromosomes, with FUL having the least and CRA having the 
most. The chromosomal locations of the DMR epimutations 
for the different finch species are shown in figure 5. All chro- 
mosomes were found to have epimutations, with CRA having 
the highest number. These chromosomal plots suggested that 
some of the species might have clusters of CNV and/or DMR 
on some of the chromosomes (figs. 3 and 4). Therefore, a 
cluster analysis previously described (Skinner et al. 2012) 
was used to examine 50-kb regions throughout the genome 
to test for statistically significant (P < 1 CT 5 ) overrepresentation 
of CNV or DMR (fig. 6). Clusters, which have an average size 
of 3 Mb, are shown as species-specific boxes for CNV (fig. 64) 
and for DMR (fig. 66). Cluster characteristics and overlap are 
presented in supplementary table S3, Supplementary Material 
online. Clusters were obtained for all species, with a higher 
number of DMR clusters than CNV clusters. The highest 
number of CNV clusters was in SCA, with more than a 4- 
fold increase over CRA (fig. 6). Therefore, in addition to 
having more CNV than expected (assuming an increasing 
number with phylogenetic distance), SCA showed more 
CNV clusters than other species (fig. 2). Genome instability 
in these cluster regions may influence the increased numbers 
of CNV in SCA, which increases the presence of CNV clusters. 
In contrast, SCA did not show more DMR numbers or clusters 
than expected, assuming an increasing number with phyloge- 
netic distance. Epimutation cluster overlap was more common 
among species (fig. 6 and table 1), suggesting that specific 
regions of the chromosomes were more susceptible to epige- 
netic alterations. Altered DNA methylation states have been 
experimentally shown to be stable for hundreds of 



Table 1 



Cluster Overlap between Species 



CNVs 


CNV 




FUL 


SCA 


PAR 


CRA 


FUL 


4 


0 


0 


2 


SCA 


0 


25 


0 


0 


PAR 


0 


0 


2 


0 


CFiA 


2 


0 


0 


6 


Epimutations 


DMR 




FUL 


SCA 


PAR 


CRA 


FUL 


16 


5 


6 


7 


SCA 


5 


16 


8 


11 


PAR 


6 


8 


16 


11 


CRA 


7 


11 


11 


25 



Note. — The overlap of CNV or DMR clusters between species is presented for 
the CNVs and epimutations. 



generations (Cubas et al. 1999; Akimoto et al. 2007; 
Skinner et al. 2010). 

The potential overlaps in specific CNV or DMR sites among 
species were examined. The overlap in genetic mutations 
among the four species is shown in a Venn diagram in 
figure 2C, whereas the overlap in epimutations is shown in 
figure 26. No overlap in specific CNV or DMR sites was ob- 
served among all species, and less than 10% overlap was 
generally observed between any two species. Interestingly, 
the CNV overlap between FUL and CRA was higher than for 
the other species (fig. 2Q. Generally, genetic and epigenetic 
alterations were distinct between species, with the majority 
being species specific. The epimutations showed more overlap 
between species than the genetic CNV mutations (fig. 26 and 
table 1). In considering within species overlap between the 
CNV and epimutations, less than 3% had common genomic 
locations. Therefore, the epimutations do not appear to be 
linked to the genetic CNV mutations, but are distinct. 

The final analysis examined the potential functional signif- 
icance of the epimutations by examining DMR and genes 
known to be associated with avian evolution. Several gene 
families and cellular signaling pathways have previously been 
shown to be involved in bird evolution, including the bone 
morphogenic protein (BMP) family and pathway (Abzhanov 
et al. 2004; Badyaev et al. 2008), the toll receptor family and 
signaling pathway (Alcaide and Edwards 201 1), and the mel- 
anins family and pathway (Mundy 2005). All the genes asso- 
ciated with these signaling pathways were localized on the 
finch genome and compared with the genomic locations of 
the epimutations and CNV. Epimutation-associated genes 
within the BMP pathway (fig. 7), toll pathway (fig. 8), and 
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Fig. 7. — Epimutation-associated genes and correlated BMP pathway. The genes having associated epimutations in the signaling pathway presented for 
the different species are identified as FUL (purple), SCA (green), PAR (blue), and CRA (red) colored boxed genes. 



melanin's pathway (fig. 9) are shown. Epimutations 
were overrepresented in all of these pathways (Fisher's exact 
test: BMP/TGFbeta (transforming growth factor) pathway, 
P< 1 x 1CT 6 ; toll pathway, P<5.7x 1CT 4 ; melanogenesis 
pathway, P< 2.5 x 1CT 13 ). Interestingly, the BMP pathway 
involved in beak development and shape had a statistically 
significant overrepresentation of CRA-associated epimutations 



when examined independently (P < 2.7 x 1 CT 5 ) (fig. 7). In ad- 
dition, the toll receptor pathway involved in immune response 
had a statistically significant overrepresentation of PAR-associ- 
ated epimutations when examined independently 
(P < 7.7 x 1 CT 4 ) (fig. 8). The melanogenesis pathway involved 
in color had a mixture of epimutations from most of the spe- 
cies when examined independently (P< 7 x 1CT 5 ) (fig. 9). 



1984 Genome Biol. Evol. 6(8): 1972-1 989. doi:10.1093/gbe/evu158 Advance Access publication July 24, 2014 



Epigenetics and the Evolution of Darwin's Finches 



GBE 




I 
& 



P P 1 

i 



5 .•; 



_i < a: < 

=> O < DC 

LL CO Q_ O 

I I I I 



Fig. 8. — Epimutation-associated genes and correlated toll receptor pathway. The genes having associated epimutations in the signaling pathway 
presented for the different species are identified as FUL (purple), SCA (green), PAR (blue), and CRA (red) colored boxed genes. 
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Fig. 9. — Epimutation-associated genes and correlated melanogenesis pathway. The genes having associated epimutations in the signaling pathway 
presented for the different species are identified as FUL (purple), SCA (green), PAR (blue), and CRA (red) colored boxed genes. 



In addition to the pathway-specific genes, the total number 
of epimutations and CNV associated with genes are presented 
in table 2, with full lists in supplementary tables S4 and S5, 
Supplementary Material online. The epimutations and CNV 
for single probe and >3 probe identification are presented 
in table 2. Observations indicate that approximately half of 
the epimutations and CNV identified were associated with 
genes. Therefore, a high percentage of the epimutations 
and CNV identified were associated with genes and were 
statistically overrepresented in several gene pathways 



previously shown to be involved in particular aspects of 
avian evolution. Although this gene association analysis dem- 
onstrates that epimutations correlate with genes and impor- 
tant pathways, the functional or causal link to specific 
evolutionary processes remains to be investigated. 

Discussion 

This study provides one of the first genome-wide comparisons 
of genetic and epigenetic mutations among related species of 
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Table 2 

Epimutation and CNV Gene Associations 



CNVs 




Total CNV 


Total CNV 


CNV Association with 


CNV Association with 




1+ Probes 


3+ Probes 


14K Genes 1+ Probes 


14K Genes 3+ Probes 


FUL 


71 


34 


40 


24 


SCA 


589 


442 


363 


350 


PAR 


295 


52 


136 


37 


CRA 


815 


602 


437 


345 


Epimutations 




Total 


Total 


Epimutation 


Epimutation Association 




Epimutations 


Epimutations 


Association with 14K 


with 14K Genes 3+ Probes 




1+ Probes 


3+ Probes 


Genes 1+ Probes 




FUL 


514 


84 


295 


48 


SCA 


890 


161 


558 


115 


PAR 


1,629 


606 


996 


407 


CRA 


2,767 


1,062 


1,611 


639 



Note. — The 14,000 zebra finch genes annotated having epimutation or CNV associations are presented for the total number of associations (overlaps) for both regions 
identified with single (1+ probes) and adjacent (3+ probes) data sets. 



organisms. There were relatively more epimutations than ge- 
netic CNV mutations among the five species of Darwin's 
finches, which suggests that epimutations are a major com- 
ponent of genome variation during evolutionary change. 
There was also a statistically significant correlation between 
the number of epigenetic differences and phylogenetic dis- 
tance between finches (figs. 1 and 3), indicating that the 
number of epigenetic changes continues to accumulate over 
long periods of evolutionary time (2-3 Myr). In contrast, there 
was no significant relationship between the number of ge- 
netic CNV changes and phylogenetic distance. 

The zebra finch genome was used as a reference for this 
study because a complete Darwin's finch genome is not yet 
available. The zebra finch genome showed hybridization with 
all probes on the array for each of the Darwin's finch species, 
suggesting that the genomes appear to be extremely similar. 
Loss of heterozygosity (absence of genomic regions, resulting 
in lack of probe hybridization) was not identified in any of the 
analyses. This suggests a high level of conservation and iden- 
tity between the species' genomes. In the event the Darwin's 
finch genome has additional DNA sequence that is not present 
in the zebra finch genome, we would not have detected this 
DNA. Therefore, our data may be an underestimate of the 
Darwin's finch genome. Another technical limitation of our 
study was that we only considered genetic CNV (amplifica- 
tions and deletions of repeat elements), but not other genetic 
variants such as point mutations or translocations. Although 
CNV frequency is higher than other mutations (e.g., SNPs) and 
stable in the genome (Gazave et al. 201 1), this study's focus 
on CNV should kept in mind. The epimutations examined are 



differential DMR that have previously been shown to be fre- 
quent and transgenerationally stable (Anway et al. 2005; 
Guerrero-Bosagna et al. 2010; Skinner et al. 2010). 
Although other epigenetic processes such as histone modifi- 
cation, altered chromatin structure, and noncoding RNA may 
also be important, DNA methylation is the most established 
heritable epigenetic mark. This aspect of the experimental 
design should be kept in mind. 

Among the five species of finches there were fewer genetic 
mutations (CNV) than epigenetic mutations. However, the 
cactus finch SCA showed a surprisingly large number of 
genetic CNV mutations than expected when compared with 
the reference species (FOR). The SCA mutations also clustered 
to similar locations on the genome to a greater extent than 
in the other species (fig. 6/4). The reason for the dispropor- 
tionately large number of CNV in the SCA comparison is 
unclear. 

In contrast to the genetic mutation (CNV) analysis, the 
number of epimutations increased monotonically with phylo- 
genetic distance (figs. 1 and 3). Overlap of specific epigenetic 
sites among species was minimal, including those for SCA (fig. 
26). An interesting possibility is that the epigenome may alter 
genome stability and generate genetic variation within spe- 
cies. A similar phenomenon has been shown for cancer, in 
which epigenetic alterations may precede genetic changes 
and alter genomic stability (Feinberg 2004). A decrease in 
the DNA methylation of specific repeat elements has previ- 
ously been shown to correlate with an increase in CNV (Macia 
et al. 2011; Tang et al. 2012). Therefore, environmentally 
induced abnormal epigenetic shifts may influence genetic 
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mutations, such that a combination of epigenetics and genet- 
ics promotes phenotypic variation. Our observations demon- 
strate a relationship between the number of epigenetic 
changes and phylogenetic distance. 

A comparison of the positions of epimutations and known 
gene families was also carried out. These gene families in- 
cluded those involved in the BMP pathway, which is related 
to beak shape (Badyaev et al. 2008), the toll receptor path- 
way, which is involved in immunological function (Alcaide and 
Edwards 2011), and the melanogenesis pathway, which af- 
fects color (Mundy 2005). Genes in all three of these families 
and signaling pathways were found to have species-specific 
epimutations (figs. 7-9). Future studies should focus on the 
causal relationship between epigenetic alterations and pheno- 
typic traits. 

Genetic mutations are postulated to provide much of the 
variation upon which natural selection acts (Gazave et al. 
2011; Stoltzfus 2012). However, genetic changes alone are 
limited in their ability to explain phenomena ranging from the 
molecular basis of disease etiology to aspects of evolution 
(Skinner et al. 2010; Day and Bonduriansky 2011; Longo 
et al. 2012; Klironomos et al. 2013). Therefore, genetic mu- 
tations may not be the only molecular factors to consider 
(Richards 2006, 2009). Indeed, epigenetic and genetic 
changes may jointly regulate genome activity and evolution, 
as recent evolutionary biology modeling suggests (Day and 
Bonduriansky 201 1; Klironomos et al. 2013). This integration 
of genetics and epigenetics may improve our understanding 
of the molecular control of many aspects of biology, including 
evolution. 

Supplementary Material 

Supplementary tables S1-S6 are available at Genome Biology 
and Evolution online (http://www.gbe. oxfordjournals.org/). 
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